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Abstract 

Canonical Polyadic Decomposition (CPD) of a third-order tensor is a minimal decomposition into a sum 
of rank-1 tensors. We find new mild deterministic conditions for the uniqueness of individual rank-1 tensors 
in CPD and present an algorithm to recover them. We call the algorithm “algebraic” because it relies only 
on standard linear algebra. It does not involve more advanced procedures than the computation of the null 
space of a matrix and eigen/singular value decomposition. Simulations indicate that the new conditions for 
uniqueness and the working assumptions for the algorithm hold for a randomly generated I x J x K tensor 
of rank R>K>J>I>2iiRis bounded as R < {I + J + K - 2)/2 + {K - P (7 - J)2 -t 4A)/2 at 
least for the dimensions that we have tested. This improves upon the famous Kruskal bound for uniqueness 
R < {I + J + K — 2)/2 as soon as 7 > 3. 

In the particular case R — K, the new bound above is equivalent to the bound R < {I — 1)(T — 1) 
which is known to be necessary and sufficient for the generic uniqueness of the CPD. An existing algebraic 
algorithm (based on simultaneous diagonalization of a set of matrices) computes the CPD under the more 
restrictive constraint 77(7? — 1) < 7(7 — 1) J(J — l)/2 (implying that 7? < (J — |)(7 — p/\/2 -|- 1). On the 
other hand, optimization-based algorithms fail to compute the CPD in a reasonable amount of time even in 
the low-dimensional case 7 = 3, J = 7, 77 = 7? = 12. By comparison, in our approach the computation takes 
less than 1 sec. We demonstrate that, at least for 7? < 24, our algorithm can recover the rank-1 tensors in 
the CPD up to 7? < (7 — 1)(J — 1). 

Keywords: canonical polyadic decomposition, CANDECOMP/PARAFAC decomposition, CP decomposi¬ 
tion, tensor, uniqueness of CPD, uni-mode uniqueness, eigenvalue decomposition, singular value decomposition 

AMS subject classifacation: 15A69, 15A23 


1 Introduction 

Let F denote the field of real or complex numbers and T € ptxjx/f denote a third-order tensor with entries 
tijk- By definition, T is rank-1 if it equals the outer product of three nonzero vectors a € F^, b € F"^, and 
c G F^: r = a 0 b 0 c, which means that tijk = cti^jCk for all values of indices. 

*Research supported by; (1) Research Council KU Leuven: Cl project cl6/15/059-nD, CoE PFV/10/002 (OPTEC), PDM 
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Dynamical systems, control and optimization, 2012-2017); (4) EU: The research leading to these results has received funding from 
the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Advanced 
Grant: BIOTENSORS (no. 339804). This paper reflects only the authors’ views and the Union is not liable for any use that may 
be made of the contained information. 
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A Polyadic Decomposition of D expresses 7" as a sum of rank-1 terms: 


R 

7 ” = 'y ' , 

r—1 


or t^Jk = y a,rb,rCkr 


r—1 


( 1 ) 


where 

a.r = [air • ■ • a/r]^ G F^, = [6ir . . . fejr]^ G F'^, Cr = [cir . . . CKr]^ € F^. 

If the number R of rank-1 terms in Q is minimal, then Q is called the Canonical Polyadic Decomposition 
(CPD) of T and R is called the rank of T (denoted by rq-). It is clear that in Q the rank-1 terms can 
be arbitrarily permuted and that vectors within the same rank-1 term can be arbitrarily scaled provided the 
overall rank-1 term remains the same. The CPD of a tensor is unique when it is only subject to these trivial 
indeterminacies. 

We write 0 as T = [A,B,C]i^, where the matrices A := [ai ... a^j;] S F^^^, B := [bi ... bj^] S 

FJxfl and C := [ci ... cr are called the first, second and third factor matrix of T, respectively. It 

may happen that the CPD of a tensor T is not unique but that nevertheless, for any two CPDs T = [A, B, 
and T = [A, B, C]j^, the factor matrices in a certain mode, say the matrices C and C, coincide up to column 
permutation and scaling. We say that the third factor matrix of T is unique. For instance, it is well known 
that if two or more columns of the third factor matrix of T have collinear vectors, then the CPD is not unique. 
Nevertheless, the third factor matrix can still be unique [71 Example 4.11]. 

The literature shows some variation in terminology. The CPD was introduced by F.L. Hitchcock in m 
and was later referred to as Canonical Decomposition (CANDECOMP) |2], Parallel Factor Model (PARAFAC) 
nans], and Topographic Components Model |21j . Uniqueness of one factor matrix is called uni-mode uniqueness 
in mi mi- Uniqueness of the CPD is often called essential uniqueness in engineering papers and specific 
identifiability in algebraic geometry papers. It is its uniqueness properties that make CPD a basic tool for 
signal separation and data analysis, with many concrete applications in telecommunication, array processing, 
machine learning, etc. I1ISII1122]. 

The contribution of this paper is twofold. First, we find very mild conditions for uniqueness of CPD and, 
second, we provide an algebraic algorithm for its computation, i.e. an algorithm that recovers the CPD from 
T by means of conventional linear algebra (basically by taking the orthogonal complement of a subspace and 
computing generalized eigenvalue decomposition (GEVD)). 

Algebraic algorithms are important from a computational point view in the following sense. In practice, the 
factor matrices of T are most often obtained as the solution of the optimization problem 


minjlf- [A,B,C]i^||, s.t. A G B G C e 


where || • || denotes a suitable norm [2d| . The limitations of this approach are not very well-known. Algebraic 
algorithms may provide a good initial guess. In Example we illustrate that even in a small-scale problem 
such as the CPD of a rank-12 tensor of dimensions 3 x 7 x 12, the optimization approach may require many 
initializations and iterations, although the solution can be computed algebraically without a problem. 

Basic notation and conventions. Throughout the paper (7* denotes the binomial coefficient, 




k\(n—k)\ ’ 

0 , 


if fc < n, 
if fc > n; 


ta, range(A), and ker(A) denote the rank, the range, and the null space of a matrix A, respectively; (the 
A:-rank of A [141 P- 162]) is the largest number such that every subset of columns of the matrix A is linearly 
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independent; “©” and denote the Khatri-Rao and Kronecker product, respectively: 


A © B = [ai 0 bi ... © bj^], 

a © b = [ai6i... aihj ... aibi... a/5j]^. 
It is well known that PD 0 can be rewritten in a matrix form as 


Ri,o(r) := 


Ti 


■BDiag(ai)C^' 

T/. 


_BDiag(a^)C^_ 


= (A©B)C^ e F 


JJxK 


( 2 ) 


where := denotes the zth horizontal slice of T = ■“ denotes the ith 

row of A £ and Diag(a®) denotes a square diagonal matrix with the elements of the vector a® on the 

main diagonal. 

To simplify the presentation and w.l.o.g. we will assume throughout that the third dimension K coincides 
with rc, yielding rc = K < R. This can always be achieved in a “dimensionality reduction” step (see, for 
instance, [3 Subsection 1.4]). 


2 Previous results, new contribution, and organization of the paper 

To explain our contribution, we first briefly recall previous results on uniqueness conditions and algebraic 
algorithms. (We refer the readers to [31313 and references therein for recent results and a detailed overview of 
early results.) 

2.1 At least two factor matrices have full column rank 

We say that a matrix has full column rank if its columns are linearly independent, implying that it cannot have 
more columns than rows. The following result is well-known and goes back to Kronecker and Weierstrass. 

Theorem 1. IISI \20^ Let T = [A, B, C]r and suppose that the matrices B and C have full column rank and 
that any two columns of A are linearly independent: 

r-B=rc= R, /ca > 2. (3) 

Then rj- = R, the CPD of T is unique and can he found algebraically. 

Theorem]^ is the heart of the algebraic algorithms presented in [3 and also in this paper. To give an idea 
of how the CPD in Theorem is computed, let us consider the particular case oi 2 x R x R tensors. Then, by 
0,B and C are Rx R nonsingular matrices. For simplicity we also assume that the second row of A does not 
contain zero entries. By ([^, PD Q can be rewritten as 

Ti = BDiag(a^)C^ and T 2 = BDiag(a^)C^, (4) 

which easily implies that 

TiT^i=BDB-\ = CT>C-\ 

where D = Diag(a^)Diag(a^) By the assumption Ica > 2, the diagonal entries of D are distinct. Hence, 
the matrices B and C can be uniquely identified up to permutation and column scaling from the eigenvalue 
decomposition of T 1 T 2 ^ and , respectively. One can then easily recover A from Q. Note that, in 

general, the matrices B and C in Theoremcan be obtained from the GEVD of the matrix pencil (Ti,T 2 ). 
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2.2 At least one factor matrix has full column rank 

In this subsection we assume that only the third factor matrix of T has full column rank. It was shown in m 
that PD 0 is unique if and only if 

rADiag(A)B^ ^ 1 for all A = (Ai,..., Xr) with at least two nonzero entries. (5) 

Condition ([^ is not easy to check for a specific tensor. The following condition is more restrictive but easy to 
check liH]. We denote by Cm(A) G the mth compound matrix of A G i.e. the matrix containing 

the determinants of all m x m submatrices of A arranged with the submatrix index sets in lexicographic order. 

Theorem 2. 0/ Let T = [A, B, C]fi and suppose that 

the matrices C 2 (A) ©C 2 (B) and C have full column rank. (6) 

Then rq- = R and the CPD of T is unique. 

It was shown in [ii that the assumptions in Theorem also imply an algebraic algorithm. The algorithm 
is based on the following relation between T and its factor matrices: 

R2,o(r) = (C2(A) © C2(B))S2,o(C)^, ( 7 ) 

in which R 2 ,o(T) denotes an CfCj x matrix whose 

((ji(2j2 - Ji - 1) - 2)/(/ - I)/4 + hi2i2 -n- l)/2, (r2 - l)i? + ri) -th 
(1 < *1 < *2 < ^, 1 < ji <j2<J, 1 < ri,r2 < R) 

entry is equal to ti.^jiriti 2 j 2 r 2 ~tiij 2 riti 2 jir 2 ~tiij 2 r 2 ^i 2 . 3 iri ^nd S 2 ,o(C) denotes an R^ x matrix 

that has columns © CrJ, 1 < ri < r 2 < R. Computationally, the identity Q is used as 

follows. First, the subspace ker(R 2 ,oCr)) is used to construct an auxiliary Rx Rx R tensor W that has CPD 
W = [C~^, C“^, M]jj in which both C~^ and M have full column rank. The CPD of W is computed as in 
Theorem[^ which gives the matrix C“^. The third factor matrix of T, C, is obtained from and the first 
two factor matrices A and B can be easily found from Ri_o(T)C“^ = A © B (see ([^) using the fact that the 
columns of A © B are vectorized rank-I matrices. 

2.3 None of the factor matrices is required to have full column rank 

The following result is known as Kruskal’s theorem. It is the most well-known result on uniqueness of the CPD. 
Theorem 3. Let T = [A,B, and suppose that 

*2R “t" 2 © T kc- (8) 

Then rp = R and the CPD of T is unique. 

In [3 H] we presented several generalizations of uniqueness Theorems and In [9] we showed that the 
CPD can be computed algebraically under a much weaker assumption than (|^. 

Theorem 4. Theorem 1.7] Let T = [A,B, and suppose that 

Cm(A) © Cm(B) has full column rank for m = R — kc + 2. (9) 

Then rj- = R, the CPD of T is unique and can he computed algebraically. 
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The algorithm in [5] is based on the following extension of Q: 

Rm,0(r) = {Cm{A) © C^{B))SmACf, (10) 

where the Cy^Cy x matrix Rm_o(T) is constructed from the given tensor T and the Cy x matrix 
Sm,o(C) depends in a certain way on C. We refer the reader to [S] for details on the algorithm. Here we just 
mention that assumption ([^ guarantees that the matrix C can be recovered from the subspace ker(Rm o(7”)). 


2.4 Generic uniqueness results from algebraic geometry 


So far we have discussed deterministic conditions, which are expressed in terms of particular A, B, C. On the 
other hand, generic conditions are expressed in terms of dimensions and rank and hold “with probability one”. 
Formally, we say that the CPD of a generic I x J x K tensor of rank R is unique if 

/r{(A, B, C) : the CPD of the tensor T = [A, B, is not unique} = 0, 
where /i denotes the Lebesgue measure on 

It is known from algebraic geometry that if2<I<J<K<R, then each of the following conditions 
implies that the CPD of a generic I x J x K tensor of rank R is unique: 

I + J + 2K -2- - jy +AK 


R < 
R < 


UK 


I + J + K-2 
IJ 


-A, 3 < /, F = C 


R < 2“+^-2 < 


(see [13 Proposition 1.6]), 

(11) 

(see [U Corollary 6.2]), 

(12) 

(see [3 Theorem 1.1]), 

(13) 


where a and /3 are maximal integers such that 2“ < / and 2^ < J. Bounds complement each other. 

If ii = A, then bound 0 is equivalent to 


R< (/-l)(J-l). 


(14) 


If F = C, then (14) is not only sufficient but also necessary, i.e., the decomposition is generically not unique for 
i? > (/ — 1)(J — 1) [3 Proposition 2.2]. 


2.5 Generic versions of deterministic uniqueness conditions 


Theorems [2||^ taken from Eiini, give deterministic conditions under which the CPD is unique and can be 
computed algebraically. Generic counterparts of condition ([^ and Kruskal’s bound ([^, for the case where 
max(/, J, A) < i?, are given by 

Cji < CjCj and R< K (see [B]) and (15) 

2i? + 2</ + J + A (trivial), (16) 

respectively. We are not aware of a generic counterpart of condition ([^, but, obviously, (|^ may hold only if 
the number of columns of the matrix Cm(A) 0 Cm(B) does not exceed the number of rows, i.e., if 


cy < cycy, where TO = i? - A + 2. 


(17) 


It can be verified that the algebraic geometry based bound (11) significantly improves bounds 0-0 if 
min(/, J) > 3. For instance, if A = A, then bound is equivalent to (14), as has been mentioned earlier, 
while ( |I^ and (16) reduce to A < (J — ^){I — \)l\/2 + 1 and A < / + J — 1, respectively. 
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2.6 Our contribution and organization of the paper 

In this paper we further extend results from PHIHIE], narrowing the gap with what is known from algebraic 
geometry. Namely, we present new deterministic conditions that guarantee that the CPD is unique and can be 
computed algebraically. Although we do not formally prove that generically the condition coincides with ( [TT] ), 
in our simulations we have been able to find the factor matrices by algebraic means up to the latter bound 
(Examples!^ and |16[ ). Moreover, the algebraic scheme is shown to outperform numerical optimization (Example 

T^. 

Key to our derivation is the following generalization of ([^, Q, and (10): 


R-m.K'T) := 4*m,i(A, B)Sm+i(C)^, m>l, Z > 0, 


(18) 


in which the matrices Rm,;(T), #m,/(A,B), and Sm+i{C) are constructed from the tensor T, the matrices A 
and B, and the matrix C, respectively. The precise definitions of these matrices are deferred to Section as 
they require additional technical notations. In order to maintain the easy flow of the text presentation, the 
proof of (18) is given in|^ The following scheme illustrates the links and shows that, to obtain our new results, 
we use (18) for to > 2 and I > 1: 


i>2.l=0 


> ( 10 ) 


n=2,Z=0 


^ (7) 


( 2 ) ( 18 ) 

j|m>2,Z>l 

Theorem 1 new results Theorem 4 => Theorem 3 Theorem 2 


(To clarify the link between (18) and (10), we need to mention that the matrices Rm^o(7”) and Cm (A) ©Cm(B) 
in ( |Io| ) are obtained by removing the zero and redundant rows of the matrices Rm_o(7”) and $m,o(A, B), 
respectively). 

Our main results on uniqueness and algebraic algorithms for CPD are formulated, explained, and illustrated 

^we generalize results mentioned 
is the special case of Theorem 15 


in Sections [4|( Theore m]^ and [^(Theorems 11-15). Namely, in Sections and 
in Subsections |2.2| and |273| respectively. In particular. Theorem]^ in Section |4 
in Section where the third factor matrix has full column rank, i.e. rc = K = R. For reasons of readability, 
in our presentation we proceed from the easy Section]^ (rc = R) to the more difficult Section [^(rc < R)- The 
proofs related to Sections and are given in Section and In Section we go from complicated to easy, 
i.e., in Subsections |6 . lH6?5] we first prove the results related to Section and then we derive Theorem from 


Theorem ]_5 in Subsection |6.6[ The paper is concluded in Section 


Our presentation is in terms of real-valued tensors and real-valued factor matrices for notational convenience. 
Complex variants are easily obtained by taking into account complex conjugations. 


3 Construction of the matrices Rto,z(T), 

Let us first introduce some additional notation. Throughout the paper denotes the set of all per¬ 

mutations of the set {li,... ,1^}- We follow the convention that if some of the values li,... ,lk coincide, then 
the cardinality of is counted taken into account multiplicities, so that always card = kl. 

For instance, ^ consists of six identical entries (1,1,1). One can easily check that any integer from 
{1,..., J"i+*} can be uniquely represented as {i — 1) +j and that any integer from {1,..., can 
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be uniquely represented as k, where 


m+/ 

P^l 

771 +/ 

+=H-^(jp-l)J™+'+ 

p^l 
771 + / 

k := 1 + ^(fcp- l)X'"+'-P, 

P^l 


*li ■ ■ • I im+l G {Ij • • ■ ; d}, 

Jl) • ■ • ) jm + l G {!) ■ ■ ■ ; d}, 
fcl, . . . , km+l G {1, . ■ . , 


(19) 

( 20 ) 
( 21 ) 


These expressions are useful for switching between tensor, matrix and vector representations. We can now 
define Il,n,i{T) as follows. 

Definition 5. Let T G . The j™+b 5 y_^"i+* matrix whose {{i — 1)7™+* + j,k)th entry is 


1 

m\(m + l)\ 


■^{^ll - +T7l + i} 




^imjmSm 


I 

n ^-^■m+pjm.+pS. 
p^l 


■+P 


( 22 ) 


is denoted by Rm^/(T). 

The matrices 4>m,i(A, B) and Sm+i(C) will have M{m,l,R) columns, where 

M(m, I, R) := C'llCZ+l, + + • • • + C™+*C™+*li. 

The columns of these matrices are indexed by (m + Z)-tuples (ri,..., rm+i) such that 


1 < n < r 2 < • • • < rm+i < R and 
the set {ri,..., r^+i} contains at least m distinct elements. 


(23) 


It is easy to show that there indeed exist M(m,l,R) (m + ^)-tuples which satisfy condition (23). We follow 
the convention that the (to + /(-tuples in (23) are ordered lexicographically: the (to -|- /(-tuple (r(,..., +_|_/) is 
preceding the (to -|- /(-tuple (r",..., + ,;( if and only if either r[ < r'( or there exists /cG{1,...,to-|-/ — 1} 


such that r( = r",... r(, = r'l and 


fe-i-i 


< r‘ 


J! 
k+l- 


We can now define $m^;(A,B( and Sm+i{C) as follows. 


Definition 6. Let A £ B £ The -by-M{m,l^R) matrix whose {{i — 1(J™+* -|- 

j, (ri,.. .,rm+i))th entry is 


Y "i®* 

^iiSi ' •■ 

^il Sm 

■ det 



■,Sm + l)6 

^imSi ■ • ■ 

^imSm 


PjmSl 





+ iSm+i 

^tm + lSm + 1 ■ 

■ ■ 


is denoted by $m,;(A,B(. 
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Definition 7. Let C G The K"^~^’‘-by-M{m,l,R) matrix whose 

in,---, rm+i)th column is 

c.„ (g) • • • « c.. 


(m + /)! ^ 

(si,....s„+,)G . 


is denoted by Sm+i{C). 


(24) 


4 At least one factor matrix of T has full column rank 


In this section we generalize results from Subsection |2.2[ i.e. we assume that the matrix C has full column rank 
and without loss of generality rc = K = R. The more general case rc = K < R is handled in Theorem in 
Section The goal of this section is to explain why and how the algebraic algorithm works in the relatively 
easy but important case rc = K = R, so that in turn Section [5] will be more accessible. 

It can be shown that for I = 0, condition (251 in Theorem l^below reduces to condition (1^. Thus, Theorem 


[^is the special case of Theorem corresponding to Z = 0. The simulations in Example below indicate that 
it is always possible to find some Z > 0 so that (251 also covers ©• Although there is no general proof, this 
suggests that ([^ can always be verified by checking (25) for some Z > 0. This would imply that Algorithm 
can compute the CPD of a generic tensor up to the necessary condition R < (I—I)(J— 1). Exampleconfirms 
this up to i? < 24. 

Let C denote the subspace spanned by all vectors of the form x (g) • • • (g) x, where 

X G is repeated m + Z times. In other words, ) contains vectorized versions of all AT x • • • x AT 

symmetric tensors of order m + Z, yielding dim We have the following result. 

Theorem 8. Let T = = [A, B,C]^, rc = K = R, I > 0, and let the matrix R 2 ,i(T) be defined as 

in Definition^ Assume that 


dim (^ker(R 2 ,;(T)) Q S 


2+lfmK^+‘- 


))=R. 


(25) 


Then 


(1) r'f = R and the CPD of D is unique; and 

(2) the CPD ofT can be found algebraically. 


Condition (25) in Theorem |8| m eans that the intersection of ker(R 2 ,i(T)) and ) has the minimal 

possible dimension. Indeed, by ( |18[ ), Definition]^ and the assumption rc = K = R, we have that the intersection 
contains at least R linearly independent vectors: 


ker(R2,,(r))f|^'+'(®'^'^') = ker(#2./(A,B)S2+i(C)^)f]52+'(M^'^') D 

ker(S 2 +/(C)^) Pi ^ ) 9 X (g) • • • (g) X, x is a column of C“^. 


The procedure that constitutes the proof of Theorem [^2) is summarized as Algorithm]^ Let us comment 
on the different steps. From Definition it follows that the rows of the matrix R 2 ,i(T) are vectorized versions 
of A" X • • • X AT symmetric tensors of order 2 + Z. Consistently, in step 2, we find the vectors -wi,... ,wji 
that form a basis of the orthogonal complement to range(R 2 ,/(T)^) in the space (K^^^*). In other words. 













spanjwi,..., Wfj} = ker(R 2 ,i(T)) nIf this subspace has minimal dimension, then its structure 
provides a key to the estimation of C. Indeed, we have already explained that the minimal subspace is given by 


ker(R 2 j(r)) f] = range 



(26) 


In steps 4-5 we recover C ^ 
matrix M such that 


from W using ([2^ as follows. By ([2^, there exists a unique nonsingular R x R 


W = (C-^ © • • • © C-^) M^. 


(27) 


In step 4, we construct the tensor W whose vectorized frontal slices are the vectors Wi,..., w/j. Reshaping 
both sides of (27) we obtain the CPD W = © ••• © In step 5, we find the CPD by 

means of a GEVD using the fact that all factor matrices of W have full column rank, i.e., we have reduced the 


problem to a situation that is covered by the basic Theorem Finally, in step 6 we recover A and B from 
R-i,o(7~)C“^ = A © B using the fact that the columns of A © B are vectorized rank-1 matrices. 


Algorithm 1 (Computation of CPD, K 


R (see Theorem 8'ii))) 


Input: T € and I > 0 with the property that there exist A S B G 

that R>2,T= [A, B, rc = R, and ( [^ holds. 

Output: Matrices A G B G and C G such that T = [A, B, C]ij 

1: Construct the / 2 +* j 2 -i-/ ^ /j 2 -i-z ^latrix R 2 .z(T) by Definition]^ 

2: Find wi,... that form a basis of ker(R 2 ,i(T)) H 
3: W •(— ["Wi . . . Wij] 

4: Reshape the x R matrix W into an R x R^+* x R tensor W 
5: Compute the CPD 

W = [C-^, C-^ © • • ■ © C-'^, M]ii (M is a by-product) (CEVD) 

6: Find the columns of A and B from the equation A © B = Ri o(T)C“^ 


and C G such 


The following example demonstrates that the CPD can effectively be computed by Algorithm for R < 
min((/ — 1)( J — 1), 24). 

Example 9. ITe consider I x J x (/— 1)(J— 1) tensors generated as a sum of R = (/— 1)(J— 1) random rank-1 
tensors. More precisely, the tensors are generated by a PD [A, B, C]i^ in which the entries of A, B, and C are 
independently drawn from the standard normal distribution N{0, 1). We try different values I = 0,1,..., until 


eondition (25) is met (assuming that this will be the case for some I > 0). We test all cases / x J x (/ — 1)( J — 1) 
such that / > 3, J > 3, and {I — 1){J — 1) < 24. The results are shown in Ta&Zej^ In all cases (25) indeed 


holds for some I < 2; the actual value of I does not depend on the random trial, i.e., it is constant for tensors of 
the same dimensions and rank. By comparison, the algebraic algorithm from mm is limited to the cases where 
I — 0, which is not always sufficient to reach the bound R < (/ — 1)(J — 1). In our implementation, we retrieved 
the vectors wi, ... ,'wr from the R-dimensional null space of a x positive semi-definite matrix 

Q. The storage of Cl is the main bottleneck in our implementation. To give some insight in the complexity of 
the algorithm we ineluded the computational time (averaged over 100 random tensors) and the size of Cl in the 
table. We implemented Algorithm^in MATLAB 2014a (the implementation was not optimized), and we did 
experiments on a computer with Intel Core 2 Quad CPUQ9650 3.00 GHzxf and 8GB memory running Ubuntu 
12 . 04.5 LTS. 
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Table 1: Values of parameter I in Theoremj^and computational cost of Algorithmfor / x J x (/ — 1)(J — 1) 
tensors of rank R = (I — 1)(J — 1) < 24 (see Example]^ for details). Note that the CPD is not generically 
unique if i? > (/ — 1)(J — 1) (see Subsection |2.4[ ). In all cases a value of I is found such that Algorithmcan 
be used. The rows with I > 1 are new results. 


dimensions of T 

/ 


computational time (sec) 

3x3x4 

0 

10 

0.02 

3x4x6 

0 

21 

0.035 

3x5x8 

0 

36 

0.051 

3 X 6 X 10 

0 

55 

0.074 

3 X 7 X 12 

1 

364 

0.403 

3 X 8 X 14 

1 

560 

0.796 

3 X 9 X 16 

1 

816 

1.498 

3 X 10 X 18 

1 

1140 

2.617 

3 X 11 X 20 

1 

1540 

5.032 

3 X 12 X 22 

1 

2024 

7.089 

3 X 13 X 24 

1 

2600 

11.084 

4x4x9 

0 

45 

0.06 

4 X 5 X 12 

1 

364 

0.401 

4 X 6 X 15 

1 

680 

1.096 

4 X 7 X 18 

2 

5985 

30.941 

4 X 8 X 21 

2 

10626 

93.03 

4 X 9 X 24 

2 

17550 

360.279 

5 X 5 X 16 

1 

816 

1.473 

5 X 6 X 20 

2 

8855 

64.116 

5 X 7 X 24 

2 

17550 

351.968 


The next example illustrates that Algorithm may outperform optimization algorithms. 
Example 10. LetT= [A,B,C]i 2 G M 3 x 7 x 12 ^ 

A = hankel{{l, 2, 3), (3, 5, 7,0, 6 , 6 , 7,9,0, 8 , 2,1)^), 

B = [I 7 hankel{{l, 2, 3,4, 5, 6 , 7), (7, 0,1, 2, 3)"^)], C = I 12 , 


where hankel{c, r^) denotes a Hankel matrix whose first column is c and whose last row is . It turns out that 
(25) holds for I = 1. It takes less than 1 second to compute the CPD ofP by Algorithm^ On the other hand, it 
proves to be very difficult to find the CPD by means of numerical optimization. Among other optimization-based 
algorithms we tested the Gauss-Newton dogleg trust region method The algorithm was restarted 500 times 
from various random initial positions. In only 4 cases the residual 

\\T — [Aest, Bes 4 , Cest]i 2 ||/||T|| after 10000 iterations was of the order 0 /0.0001 and in all cases the estimated 
factor matrices were far from the true matrices. Other optimization-based algorithms did not yield better results. 
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5 None of the factor matrices is required to have full column rank 

In this subsection we consider the PD T = {tijkYij'k=i — and extend results of the previous 

subsection to the case rc = K < R. 


5.1 Results on uniqueness of one factor matrix and overall CPD 

We have two results on uniqueness of the third factor matrix. 

Theorem 11. Let T = = [A, B, C]/j, rc = K<R, m = R — K + 2, and li,... ,lm be nonnegative 

integers. Let also the matrices (A, B),..., (A, B) and Si_|_/j^(C), Sm+i„(C) be defined as in 

Definition and Definition respectively. Let Ui,..., be matrices such that their columns form bases for 
range(Si+q (C)^),..., range(Sm+i„ (C)^), respeetively. Assume that 

(i) fee > 1; and 

(ii) A © B has full column rank; and 


(hi) (A, B)Ui, (A, B)Um have full column rank. 


Then rp = R and the third factor matrix ofT is unique. 


According to the following theorem the set of matrices in (iii) in Theorem 
matrix if i? < min(fcA, ks) + K — 1. 


11 can be reduced to a single 


Theorem 12. Let T = {tijkYfif'^i ~ [A, B,C]i^, rc=K<R, m = R—K + 2, and I > 0. Let also the 
matrices (A,B) and Sm+i(C) be defined as in Definition^ and Definition^ respeetively. Let Um be a 
matrix such that its columns form a basis for ia.nge{Sm+i{C)'^). Assume that 


(i) kc > 1; and 


(ii) A 0 B has full column rank; and 


(iii) min(/cA) ^b) > w — 1; and 


(iv) the matrix $m^;(A, B)Um has full column rank. 
Then rp = R and the third factor matrix ofT is unique. 


The assumptions in Theorems and complement each other as follows: in Theorem 11 we do not 
require that the condition min(fcA) ^b) > m — 1 holds while in Theorem 12 we do not require that the matrices 
^k,ik (A, B)Ufe, 1 < fc < TO — 1 have full column rank. 

It was shown in [SI Proposition 1.20] that if T has two PDs T =[A, B, Cj^j and T = [A, B, C]fl that share 
the factor matrix C and if the condition 


max(min(fcA, k^ — 1), min(fcA — 1, ^b)) + kc > R + 1 


(28) 


holds, then both PDs consist of the same rank-one terms. Thus, combining Theorems 11 -12 with |8l Proposition 
1.20] we directly obtain the following result on uniqueness of the overall CPD. 


Theorem 13. Let the assumptions in Theorem 11 or Theorem 12 hold and let condition (28) be satisfied. Then 
rp = R and the CPD of tensor T is unique. 
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5.2 Algebraic algorithm for CPD 

We have the following result on algebraic computation. 

Theorem 14. Let T = ~ rc=K<R, m = R—K +2, and I > 0. Let also the 

matrices ^m,i (A,B) and Sm+i(C) be defined as in Definition^ and Definition^ respectively. Let Um be a 
matrix such that its columns form a basis for range(Sm+;(C)^). Assume that 

(i) kc = K; and 

(ii) A © B has full column rank; and 


(ill) the matrix B)Um has full column rank. 

Then rj- = R, the CPD of T is unique and can be found algebraically. 

The assumptions in Theorem are more restrictive than the assumptions in Theorem as will be clear 
from Section Hence, the statement on rank and uniqueness in Theorem follows from Theorem 13 To 
prove the statement on algebraic computation we will explain in Section that Theorem 14 can be reformulated 
as follows (see Section]^ for the definition of ^ )). 

Theorem 15. LetT = {U 3 k)['j’k=i = [A,B,C]_r, rc = K < R, m = R — K-\-2, and 1 > 0. Let also the matrix 
be defined as in Definition^ Assume that 

(i) kc = K; and 

(ii) A © B has full column rank; and 


(hi) dim (^ker(R™,z(T)) fl S 




))=C 


K-1 


Then rj- = R, the CPD of T is unique and can be found algebraically. 


Note that if kc = K, then by (18) and Lemma |2^ (i) below, 


dim (ker($™AA, B)S^+z(C)^) f] 
dim (ker(S™+z(C)^) = C 


> 


(29) 


K-l 

fl 


Thus, assumption (iii) of Theorem 15 means that we require the subspace to have the minimal possible dimen¬ 
sion. That is, we suppose that the factor matrices A, B, and C are such that the multiplication by $m,i(A,B) 
in (18) does not increase the overlap between ker(S m+iiC)'^) and ^ ). In other words, we suppose 

that the multiplication by 4>m,i(A, B) does not cause additional vectorized K x ■ ■ ■ x K symmetric tensors of 
order m + I to he part of the null space of Rm,i(T). This is key to the derivation. By the assumption, as we 
will explain further in this section, the only vectorized symmetric tensors in the null space of Rrn.iiT) admit a 
direct connection with the factor matrix C, from which C may be retrieved. On the other hand, the null space 
of R„i_/(T) can obviously be computed from the given tensor T. 

The algebraic procedure based on Theorem consists of three phases and is summarized in Algorithm 
In the first phase we find the K x C^~^ matrix F such that 


every column of F is orthogonal to exactly K — 1 columns of C and 


(30) 
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any vector that is orthogonal to exactly K — \ columns of C 

is proportional to a column of F. 


( 31 ) 


Since kc = K any K — 1 columns of C define a unique column of F (up to scaling). Thus, (30)-(31| define 


Algorithm 2 (Computation of CPD, K < R (see Theorem 

15 

1) 

Input: T and 1 > 0 with the property that there 63 

that T = [A, B,C]i? and assumptions (i)-(iii) in Theorem 
Output: Matrices A G B G and C G such 

dst A G B G and C G R^^^ such 

15 hold. 

that T = [A, B, Cji,; 

Phase 1: Find the matrix F G such that F coincides with B{C) up to (unknown) column permutation 

and scaling 


1: Construct the jm+ijm+i X j^m+i jnatrix R™,z(r) by Definition^ 

2: Find Wi,..., WpK-i that form a basis of ker(Rm^;(T)) H •S'™'*'*) 

3 : W ^ [wi . . . 

4: Reshape the x matrix W into an RT x x tensor W 

5: Compute the CPD 

W = [F, F © • • ■ © F, (M is a by-product) (GEVD) 

Phase 2 and Phase 3 (can be taken verbatim from [9l Algorithms 1,2]) 


the matrix F up to column permutation and scaling. A special representation of F (called B{C)) was studied 
in [9]. It was shown in [9] that the matrix F can be considered as an unconventional variant of the inverse of C: 

every column of C is orthogonal to exactly C^Zi columns of F, 
any vector that is orthogonal to exactly C^Zi columns of F 

is proportional to a column of C. 


(32) 

(33) 


(Note that, since fee = K, multiplication by the Moore-Penrose pseudo-inverse yields CC^ = T-k- In 
contrast, for F we consider the product FC.) It can be shown (see Lemma 231 that under the assumptions in 
Theorems [MKIsI 


kjp > 2 , 


the matrix has full column rank and 


where 


ker(R™,,(r)) f| = range (f('"+')) , 


p(m+i-i) .= F © • • • © F, := F © • • • © F . 


(34) 

(35) 

(36) 


1 


m+Z 


If A = i? (as in Subsection, then m = R—K + 2 = 2, (35) coincides with (26) (F coincides with C ^ up 


to column permutation and scaling), and the first phase of Algorithm ^coincides with steps 1-5 of Algorithm 
For K < R (implying m > 2) we work as follows. From Definition it follows that the rows of the matrix 
R-m,i(T) are vectorized versions oi K x ■ ■ ■ x K symmetric tensors of order m 1. Thus, in step 2, we find 
the vectors Wi,... that form a basis of the orthogonal complement to range(R,„_i(T)^) in the space 
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^ existence of such a basis follows from assumption (hi) of Theorem 15). By (35), there exists 
a unique nonsingular 


K-l 


X ^ matrix M such that 


■Yv = 


(37) 


In step 4, we construct the tensor W whose vectorized frontal slices are the vectors Wi,... ,w^k-i. Reshaping 

both sides of (37 ) we obtain the CPD W = [F, in which the matrices and M have full 

column rank and kp > 2. By Theorem the CPD of W can be computed by means of GEVD. 

In the second and third phase we use F to find A, B, C. There are two ways to do this. The first way is to 
find C from F by ([M|)-([M| and then to recover A and B from T and C. The second way is to find A and B 
from T and F and then to recover C. The second and third phase were thoroughly discussed in and can be 
taken verbatim from [9l Algorithms 1 and 2]. 

Example 16. Table contains some examples of CPDs which can be computed by Algorithm and cannot 
be computed by algorithms from m- The tensors were generated by a PD [A, B,C]ij in which the entries of 
A, B, and C are independently drawn from the standard normal distribution A(0,1). We try different values 
I = 0,1,..., until condition (hi) in Theorem 15 is met (assuming that this will be the case for some I > 0). 
In our implementation, we retrieved the vectors wi,..., from the -dimensional null space of a 

xC'fl+/+i positive semi-definite matrix Q. The storage o/Q is the main bottleneck in our implementation. 
To give some insight in the complexity of the algorithm we included the computational time (averaged over 100 
random tensors) and the size of Q in the table. 


Uniqueness of the CPDs follows from Theorem 15 By comparison, the results of m guarantee uniqueness 
only for rows 1-4 (see Table 3.1]). 


Table 2: Upper bounds on R under which the CPD of a generic Ix JxK tensor can be computed by Algorithm 
(see Example [T^ details). 


dimensions of T 

R 

m 

1 


computational time (sec) 

4x5x6 

7 

3 

1 

126 

0.182 

5x7x7 

9 

4 

1 

462 

1.598 

6x9x8 

11 

5 

1 

1716 

28.616 

7x7x7 

10 

5 

1 

924 

8.192 

4x6x8 

9 

3 

1 

330 

0.63 

4 X 7 X 10 

11 

3 

1 

715 

2.352 

5x6x6 

8 

4 

2 

462 

1.256 

5x7x8 

10 

4 

2 

1716 

14.552 


6 Proofs related to Sections [4] and 


In this section we 1) prove Theorems [ll| and [I^ 2) show that the assumptions in Theorem 14 


are more restrictive 


than the assumptions in Theorem 13 which implies the statement on uniqueness in Theorem 14 3) prove that 
assumption (iii) in Theorem 14 is equivalent to assumption (iii) in Theorem |15[ 4) prove statements (34)-(35); 
5) prove Theorem]^ 
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6.1 Proofs of Theorems [TT] and [T^ 

In the sequel, a;(Ai,..., Xr) denotes the number of nonzero entries of [Ai ... Aj^]^. The following condition 
(Wm) was introduced in [71IH] in terms of m-th compound matrices. In this paper we will use the following 
(equivalent) definition of (Wm) . 

Definition 17. We say that condition (Wm) holds for the triplet of matrices (A, B, C) S x x 

if uj{Xi ,..., Xr) < m — 1 whenever 

rADiag(Ai.....AH)Br’< TO-1 and [Ai ... S range(C'^). (38) 

Since the rank of the product ADiag(Ai,..., A_r)B^ does not exceed the rank of the factors and rDiag(Ai,...,AH) = 
a;(Ai,..., Xr), we always have the implication 


j{Xi,...,Xr) < m - 1 


rADiag(Ai,...,AH)Br’ < m — I. 


(39) 


By Definition 17, condition (Wm) holds for the triplet (A, B, C) if and only if the opposite implication in (39) 
holds for all [Ai ... A_r] S range(C’^) C R^. 

The following results on rank and uniqueness of one factor matrix have been obtained in [7]. 

Proposition 18. (see Proposition f.Q]) Let T = = [A,B, C]i^, rc = K < R- Assume that 

(i) kc> 1; 

(ii) A © B has full column rank; 

(Hi) conditions (Wm),..., (Wi) hold for the triplet of matrices (A, B, C). 

Then rp- = R and the third factor matrix ofT is unique. 

Proposition 19. (see Corollary 4-10]) Let T = {Ujk)l'j'k=i = CJij, rc = K < R. Assume that 

(i) kc> 1; 

(ii) A © B has full column rank; 

(Hi) min(fcA, fce) > m — 1; 

(iv) condition (Wm) holds for the triplet of matrices (A,B,C). 

Then rq- = R and the third factor matrix ofT is unique. 


One can easily notice the similarity between the assumptions in Theorems |11H12| and the assumptions in 
Propositions [TsHT^ The proofs of Theorems [TT}fT^ follow from Propositions fTMT^ and the following lemma. 


Lemma 20. Let A G B G and C G rc = K < R, k<m = R — K + 2, and let I be a 

nonnegative inte^r. Let also the matrix $fc^;(A, B) be defined as in Definition^ the matrix Sk+i{C) be defined 
as in DefinitionVA and U fee a matrix such that its columns form a basis for range(Sj,+;(C)^). Assume that 


the matrix $fc,i(A, B)U has full column rank. 
Then condition (Wk) holds for the triplet of matrices (A,B,C). 


(40) 
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Proof. Let (38) hold for m = k. We need to show that Cli(Ai, ..., Xr) < k — 1. Since [Ai ... Ai^]^ G range(C^) 
and re = K, there exists a unique vector x G such that [Ai ... Ai?] = x^C. Hence, we need to show that 
X is orthogonal to at least R — k + 1 columns of C. 


By (38), there exist A G and B G such that 

ADiag(Ai,..., Ai?)B'^ = AB^ 


(41) 


and max(r^, ) < k—1. Since ab^ A = a®b®A, we can consider (41) as an equality of two PDs of an / x J x 1 

tensor 


'y ] a^ ® br ® Af = y ] a^ ® br ® 1. 


r—1 


r—1 


Hence, by (18), 


$fci(A,B)Sfc+z(x^C)^ = $fe,i(A,B)Sfe+i([Ai ... Ai?])^ = 

$fe,,(A,B)Sfe+i([l ... 1])^. 


(42) 


Since max(r^,r 3 ) < fc — 1, it follows from Dehnition that 4>fc,i(A, B) is the zero matrix (cf. explanation at 
the end of Section |^. Besides, it easily follows from Definition that 

Sfc+/(C)^(x(g)- -(g)xJ = Sfc+/(x^C)^. 
k-\-l 


Thus, (42) takes the form 

#fe,i(A, B)Sfc+,(C)^(x 0 ... 0 x) = $fc,j(A, B)Sfc+j(x^C)^ = 0. 
Hence, by (40), the vector x 0 • • • (g) x is orthogonal to the range of Sk+i {C). In particular. 


E 


-S^ + l 


(si, 


■ Sk + l)^ 


(x^c^i) • • • (x^c^^_ j(x^c^j'+i = 0 

for all {k + Z)-tuples (ri,..., r^,..., r^) such that 1 < ri < • • • < r/j < i?, yielding that x is orthogonal to at 
least i? — fc + 1 columns of C. □ 


6.2 Proof of statement on rank and uniqueness in Theorem 14 


In Lemma 21 below we prove that min(fcA) ^b) > w. It is clear that condition min(fcA, ^b) > w and assumption 
(i) in Theorem 14 imply assumption (hi) in Theorem 12 and condition (28). Hence, by Theorem 13 rq- = R 
and the CPD of tensor T is unique. 


Lemma 21. Let assumptions (i) and (hi) in Theorem 14 hold. Then 
min(fcA, ^b) > m. 
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Proof. Assume to the contrary that < m or k-Q < m. W.l.o.g. we assume that the first m columns of 
A are linearly dependent. We will get a contradiction with assumption (iii) by constructing a nonzero vector 
f G range(Sm+z(C)^) such that B)f = 0. Since kc = K, there exists x G such that 

X^Ci 7^ 0, 7^ 0, x'^C^+i = • • • = x'^Cij = 0. (43) 


We set f = S„+,(C)^(> 


c) and we index the entries of f by (m + Z)-tuples as in (23). One can easily 


m-\-l 


show that f has entries ... ( 


X Cr 


,). Hence, by (43), 


(x^c^J ... (x^Cr„+,) = 0, if {ri,... ,r™+i} \ {!,.. .,m} 7 ^ 0, 

(x'^c^J ... (x'^Cr^^,) 7 ^ 0, if {ri,... ,rr„+i} \ {1,... ,m} = 0. 

On the other hand, by Definition and the assumption of linear dependence of the vectors ai,...,am, the 
columns of $m_;(A,B) indexed by the {m + Z)-tuples (23) such that {ri,... ,rm+i} \ {1,... ,m} = 0 are zero. 
Hence, B)f = 0. □ 


6.3 Properties of the matrix Sm+z(C)" 


The following auxiliary Lemma will be used in Subsections 6.4 and 6.5 Since the proof is rather long and 
technical, it is included in[Bl 


Lemma 22. Let C G kc=K,m = R — K + 2,l>0, letF satisfy (301 -(31), and let be defined 

by (36). Then 

(i) dim(ker(S„+i(C)^)n^’”+'(K^™^')) = 

(ii) ker(S„+i(C)^)n^'"+'(K^'"^') = range (f(™+')); 


(iii) range(S„+i(C)^) = S^+,(C)^(^™+'(K^"‘^')); 

(iv) dim (range(S™+i(C)^)) = . 


6.4 Proof of equivalence of Theorems 14 and 15 


We prove that assumption (iii) in Theorem is equivalent to assumption (iii) 
sufficient to prove that 


in Theorem 15 


By (29), it is 


dim fker(R™,i(T)) Q S' 


m+l/mK^+‘ 


the matrix $m,i(A,B)U„ 




1 


does not have full column rank. 


(44) 


To prove (44) we will use the following result for X := $m,z(A,B), Y := S^+,(C)^, and E := ): 

if if is a subspace and X and Y are matrices such that XY is defined, then 


dim (ker(XY) n if) > dim (ker(Y) n if) + 1 

there exists a nonzero vector f G if \ ker(Y) such that XYf = 0. 


(45) 
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We have 




>Cn 


K-1 


1 = 


dim (ker(R„,,(r)) 

dim (ker($„,i(A, B)S„+i(C)^) f| 

dim (ker(S™+,(C)^)f|^™+'(K^’"^')) + 1 

j there exists a nonzero vector f e ^ ) \ ker(Sm+i(C)^) 1 

\ such that $m,i(A,B)Sm+;(C)'^f = Oj 

the matrix $m^;(A, B)Um does not have full column rank. 


where the equality in the second statement holds by Lemma 22 (i) and the last equivalence follows from 
range(Um) = range(S„+;(C)^). 


6.5 Proof of the statement on algebraic computation in Theorem 14 


The overall procedure that constitutes the proof of the statement on algebraic computation is summarized in 
Algorithm and explained in Subsection 5.2 In this subsection we prove statements (34)-(35|. 


Lemma 23. Let assumptions (i) and (iii) in Theorem 15 hold and let F satisfy (30)-(31). Then (34)-(35) 
hold. 


Proof. The implication kc = K ^ (34) was proved in 0 Proposition 1.10]. In Subsection 6.4 we proved that 
assumption (iii) in Theorem 14 holds. By ([Ts]), Theorem [l4| (iii), and Lemma 22 we have 


ker(R™,(r))n5'"+'(R'^™^') = 

kevi^m.i (A, B)S^+, (C)^) f| ) = 

keiiSm+iiCf ) n 5™+'(K^"'^') = range 


which completes the proof of (35). 


□ 


6.6 Proof of Theorem! 


We check the assumptions in Theorem 15 for m = 2. Assumption (i) holds since rc = K = R implies 


kc = K and assumption (iii) coincides with (25). To prove assumption (ii) we assume to the contrary that 


(A©B)[Ai ... — 0. Then rADiag(Ai,...^)Br) — 0. In Subsection 6.4 we explained that assumption (iii) 


in Theorem 14 also holds. Hence, by Lemma 20 condition (W2) holds for the triplet (A, B, C). Hence, at most 


one of the values Ai,..., A^ is not zero. If such a Xr exists, then = 0 or = 0 yielding that min(fcA, ^b) =0. 
On the other hand, by Lemma 21 min(fcA) ^b) > 2, which is a contradiction. Hence, Ai = • • • = Xr = 0. 


7 Discussion 

A number of conditions (called (Km), (C m ), (u m ), and (Wm)) for uniqueness of CPD of a specific tensor have 
been proposed in [Hi. It was shown that each subsequent condition in (Km),..., (Wm) is more general than 


18 























the preceding one, but harder to use. Verification of conditions (Km) and (Cm) reduces to the computation 
of matrix rank. In contrast, conditions (Um) and (Wm) are not easy to check for a specific tensor but hold 
automatically for generic tensors of certain dimensions and rank |10j . 

In this paper we have proposed new sufficient conditions for uniqueness that can be verihed by the com¬ 
putation of matrix rank, are more relaxed than (Km) and (Cm), but that cannot be more relaxed than (Wm). 
Nevertheless, examples illustrate that in many cases the new conditions may be considered as an “easy to check 
analogue” of (U2) (<^ (W2)) and (Wm). 

We have also proposed an algorithm to compute the factor matrices. The algorithm relies only on standard 
linear algebra, and has as input the tensor T, the tensor rank R, and a nonnegative integer parameter 1. 
The algorithm basically reduces the problem to the construction of a x matrix Q, the 

computation of its ^-dimensional null space, and the GEVD of a x matrix pencil, where 

m = R — K + 2. For I = 0, Algorithms and coincide with algorithms from [S] and [2, respectively. 
Our derivation is different from the derivations in [6] and [9] but has the same structure: from the CPD 
T = [A, B, C]/j we derive a set of equations that depend only on C; we find C from the new system by means 
of GEVD, and then recover A and B from T and C. 

It is interesting to note that the new algorithm (with I = 1) computes the CPD of a generic 3 x 7 x 12 tensor 
of rank 12 in less than 1 second while optimization-based algorithms (we checked a Gauss-Newton dogleg trust 
region method) fail to find the solution in a reasonable amount of time. 

We have demonstrated that our algorithm (with I < 2) can find the CPD of a generic I x J x K tensor of 
rank R\i R< K < (/ — 1)( J — 1) and R < 24. We conjecture that the algorithm (possibly with I > 3) can also 
find the CPD for R > 25. (It is known that the CPD of a generic tensor is not unique if i? > (/ — 1)( J — 1)). 
In that case the x C'^Xln+i-i matrix Q becomes large and the computation, as it is proposed in the 

paper, becomes infeasible. Since the null space of Q is just i?-dimensional the approach may possibly be scaled 
by using iterative methods to compute the null space. 


A Derivation of identity (18) 


Let T = {tijk)i’j 'kLi = [-A-,B, C]ii;. In this section we establish a link between the matrix Rm^;(T) defined 
in subsection and the factor matrices A, B, and C. We show that the matrix HmjiT) is obtained from 
T by taking the following steps: 1) taking the (m -I- /)th Kronecker power of T; 2) making two partial skew- 
symmetrizations and one partial symmetrization of the result; 3) reshaping the result into an j™+* x 
matrix. The main identity is obtained by applying steps l)-3) to the both sides of Q. 


A.l Step 1: Kronecker product power of T 

The Kronecker product square of T, := T 0 T, is an / x J x AT block-tensor whose {i,j, k)th block is the 

I X J X K tensor tijkT- Equivalently, is an P x J'^ x tensor whose {i,j, k) := {ii — l)I + i 2 , (ji — 1)=/ + 

72, (fci — 1)A -I- fc2)th entry is 

hjk ~ 

Similarly, the (Z -I- m)-th Kronecker product power of T, 


.= 7-^... ^7- 

m-\-l 
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is an /™+* X x tensor whose (i,j, fc)th entry is 


ijk 


= t 


iljlkiti2j2k2 


n + ljm+lkm + l I 


(46) 


R 


where i ,j, and k are defined in (191, (20), and (21), respectively. One can easily check that if T = X) a,. ®br®Cr 

r—1 

then 

R 


j-{m+l) = ^ (£ 


^r^+i) ® (bri 0 ■ ■ • 0 ) ® (Cri ® ■■■ ® ). 




A.2 Step 2: two partial skew-symmetrizations and one partial symmetrization of 
a reshaped version of 

Recall that a higher-order tensor is said to be symmetric (resp. skew-symmetric) with respect to a given group of 
indices or partially symmetric (resp. skew-symmetric) if its coordinates do not alter by an arbitrary permutation 
of these indices (resp. if the sign changes with every interchange of two arbitrary indices in the group). 

Let us recall the operations of (complete) symmetrization and skew-symmetrization. With a general A:th- 
order L x ■ ■ ■ x L tensor Af one can associate its symmetric part (Af) and skew-symmetric part (Af) as 

follows. By construction, is a tensor whose entry with indices h, • ■ •, h- is equal to 

— ^ ?T-pi...pfc- (47) 

That is, to get S^{AJ') we should take the average of k\ tensors obtained from Af by all possible permutations 

of the indices. Similarly, A^{AJ') is a tensor whose entry with indices li,... ,1^ is equal to 


H E sgn(pi,... ,pfc)npi..,p,, 

0 . 


if h) • ■ • dfe Eire distinct, 
otherwise, 


(48) 


where 

sgn(pi,... ,pk) denotes the signature of the permutation (pi,... ,pk). 

The definition of A^{AA) differs from that of S^{AA) in that the signatures of the permutations are taken into 
account and that the entries of A^(A/') with repeated indices are necessarily zeros. One can easily check that if 
AC = di ® ... ® dfe (that is, Af is the A:th order rank-1 tensor), then 


5'=(di ®. 

.. ® dfc) = 

1 

(pi.. 

^ dp, ® . 

■■,Pfe)6P{l,...,fc} 


(49) 

A'=(di ®. 

.. ® dfc) = 

1 

H 

(pi.. 

..,Pfe)eP{i,...,fe} 

■ ^Pk)^pi ® ® . 

(50) 


Partial (skew-)symmetrization is a (skew-)symmetrization with respect to a given group of indices. Instead of 
presenting the formal definitions we illustrate both notions for an M x L x L tensor J\f = {nmhi 2 )m'h’i 2 =i- 
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Partial symmetrization with respect to the group of indices {2,3} maps the tensor Af to a tensor that we denote 
by (Im ® whose entry with indices (m, hjh) is equal to 


(pi,P2)e-P{ii,i2} 

Similarly, by (Im ® A^)A/’ we denote the tensor whose entry with indices (m, ^1,^2) is equal to 


I] sgn(pl,p2)^^n^^l^2 = rimhh - nmhh, if h ^ h, 

(pi.P2)e-P{ii,i2} 

0, if li = ^2- 


If Af = di ® d2 ® da e rMxLxl^ 

(Im ® <S'^)(di ® da ® da) = di ® S'^(d2 ® da) = di ® da ® da + di ® da ® da, (51) 

(Im ® A^)(di ® da ® da) = di ® A^(d2 ® da) = di ® da ® da — di ® da ® da. (52) 

Thus, operations (Im ® S'^) and (Im ® A^) symmetrize and skew-symmetrize the horizontal slices of A/”. 

Let us reshape the tensor an Ix---xIxJx---xJxKx---xK (each letter is repeated 

m + I times) tensor as: 

R 

^ (a^i ® ... ® ® (bri ® ... ® ® (Cri ® ... ® Cr„_^,). (53) 


Then the entries of are given by 

■^m+V) 


n + in--- 3 m + lKl. 


~ fjlilfelf*2i2fc2 ■ ■ ■ ^im + ljm + lkm+r 


(54) 


From (46) and (54) it follows that T^l^+O is just a higher-order representation of 

A new tensor is obtained from 7^(™+b applying two partial skew-symmetrizations and one partial 

symmetrization as follows: 


,_ \ ( \ri 

'AAS ip'¬ 


ll/) ® (A™ ® Ij ® ... ® ij) ® 5'™+'] 


(55) 


To obtain we first skew-symmetrize with respect to the group of indices {!,...,m} (the 

first m dimensions), then we skew-symmetrize the result with respect to the group of indices 
1,..., 2m + 1} (the first m “J” dimensions), and, finally, we symmetrize the result with respect to the group 
of indices {2m -|- 2^ -|- l,...,3m -|- 3^} (all “iP” dimensions). From (47), (48), and (54), it follows that the 


{m+l) 


{ii,... ji,..., jm+h ki, ■ ■ ■, km+i)th entry of the tensor Tj^j^g 

in ii,..., Zm or ji,..., jVn and is equal to 


is equal to zero if some index is repeated 


1 

(m + l)\ 


E 

(si,...,s„,+i)G 


—T E sgn(9 i,...,g^)x 

ml 

L (gi,...,q^)G 
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^ ^ Sgn(pi, . . . ,Pm) 




(pi,...,Pm)6 

P{.l. 


1 


(m + ly. 


(si,...,s„+i)G 

p, 


E 

..,S^ 

1 

ml 


E E sgn(9i,...,g^); 

m! 


-^{31.im} 


det 


tiiqiSi 


1 




pimQlsi 


V — 1 


^ + 0! 


E 


det 


(si,....Sm + l)G 
P{k-i,...,k^^l} 


kiijisi 


pirnjlSi 


‘'llJmSr, 


11' 

v—1 


lm+v3m + vSm+v 5 


otherwise (we used twice the Leibniz formula for the determinant). Thus, by (22), the tensor T/Pjpg^'^ and the 
matrix Rm^;(T) have the same entries (in step 3 it will be shown that Rm_;(T) is a matrix unfolding of T'jPPpg^'^). 

Let us apply partial skew-symmetrizations and partial symmetrization to the right-hand side of (531: from 
(53), (55), (see also (51)-(52) for the properties of the outer product) it follows that 


R 


^ [A"*(arj ® ... ® a^^) ® SLr^^i ® ... ® a^ 
A"* (bj-i ® ... ® ) ® br 


^Tm + l ® ■ ■ ■ ® J ® 5 '” 


(56) 




where the expressions A'"(arj ® ... ®ar„) and A"*(brj ® ... ®br^) are defined in (50), the expression (c^.^ ® 


’ is defined in (|49[), and, by definition. 




A,B 


ri,...,r^+i :=A'"(ari ® ... ® ar„) ® a^^^^ « 
K^{hri ® . . . ® b^^) ® hr^^i 


-’r-m+1 


(recall that the vectors a^ and b^ are columns of the matrices A and B, respectively). 


Note that by construction, S'"*+*(cri ® ... < 


S^+^{Cri ® ... 


,) is a completely symmetric tensor ( that is, the expression 


,) does not change after any permutation of the vectors Cn ,..., Cr^ ^). Taking this fact 


into account we can group the summands in (56) as follows 


'AAS 


E 


E 




A,B 


Pl,---.Pm + l 


. ^m + l} 


s"^+ycr 


(57) 
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into the matrix Rm,z(T) 


A.3 Step 3: Reshaping (nnfolding) of 

We define the matricization operation 

Matr ■ 1^7’"+'j’"+'xX’"+* 


as follows: the (ii,..., ■ • • ,jm+i, fci, ■ ■ ■, entry of a tensor is mapped to the ((i — 1) fc)th 

entry of a matrix, where i, j, and k are defined in (19l, (20), and (21), respectively. One can easily verify that 

i b. 


Matr 

[an 


^Im + l 

a 


-*3i 


+l 


-'Jm + l 

5b 


® . . . ® = 




(Cfci O ■ 




(58) 


and that 'Rm,i{T) = Matr{T^^g^'’). 

What is left to show is that the matricization of the right-hand side of (57) coincides with the matrix 
Sm+i{C)'^■ In the sequel, when no confusion is possible, we will use S'^ and to denote “sym- 
metrization” and “skew-symmetrization” of vector representations of a certain tensor: if di,... ,dfe G K^, then 
the vectors S'^(di 0 • • • 0 d/j) and A^(di 0 • 

“®” replaced by “O”: 


I dfe) are computed in the same way as in (49)-(50) but with 


5''=(di0-. 

• • 0 dfe) = 

1 

H 

(pi.. 

E dp,0- 

••0dp,, 


A'=(di0-. 

• • 0 dfe) = 

1 

(pi.. 

E crbl.-- 

■■,Pfc)6f’{l,...,fc} 

• ,Pfc)dpi 0 • • • 0 dp^. 

(59) 


Hence, by (56), (57), and (58) 


R„M(r) = Matr(firs^')) = 

H ^ fn*5™+'(cri = 




(60) 


where 


fA,B 


:=A'"(a, 


(g) a^^) (g) as^^j 


a + l'^ 


A'"(bsi 0 • • • 0 b^^) 0 bs^^i 0 • • • 0 bs^^,. 




- E 


fA,B 


(si,...,s„+i)G 


We show that (/<(A, B)ri,...,r„+i is the zero vector if the set {ri,..., rm+i} has fewer than m distinct elements 
and that (/)(A, B)ri,...,r„+i is a column of 4>m,i('^iR) otherwise. From (59) and the Leibniz formula for the 
determinant it follows that the entries of the vector A^(di 0 • • • 0 d^) are all possible kx k minors of the matrix 
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D := [di ... dfc] divided by k\. In particular, if some of the vectors d.^ coincide, then A*^(di 0 • • • 0 d^) is the 
zero vector. Hence, the vector „ , has entries 

7 ,*rn + I 


1 

(m!)2 



Si • • • 

^2lSm 


bjisi 


det 

• • • 


• det 

PjmSl 









where G and G {1, In particular, if the set {ri,..., r^+i} has 

fewer than m distinct elements, then „ ,, are zero vectors for all (si,..., Sm+z) G P{ri,...,r^+i}, yield¬ 


ing that (j){A, B)r 
(^(A, 'R)r-i,...,rm+l 1 


r„+i is the zero vector. Hence, by Definition Ibl the matrix $m.z(A,B) has columns 
where (ri,... ,rm+i) satisfies (23). Thus, (60) coinci^s with (|l8|. 


B Proof of Lemma 22 


The proof of Lemma is based on the following simple generalization of the rank-nullity theorem and relies 
on two bounds that will be obtained in in |B.l| and |B.2[ respectively. 

Lemma 24. Let "K be a matrix and E be a subspace such that range(X^) C E. Then 

(i) dim(ker(X) O E) + dim(X(i?)) = dim£l; 

(ii) range(X) = X{E), 

where the subspace X(i?) denotes the image of E under X. 

Proof. Let P be a matrix whose columns form a basis for the subspace E. Then dim(ker(X)nill) = dim(ker(XP)), 
X(£l) = range(XP), and the matrix XP has dim£l columns. Hence, by the rank-nullity theorem, 


dim(ker(X) n i?) -I- dim(X(i5)) = dim(ker(XP)) -|- dim(XP) = dimi?. 
Since, range(X^) C range(P), it follows that range(X) = range(XP) = X(iil). 


□ 


Proof of Lemma 22. We set X = Sm+z(C)^ and E = ). Then statement (iii) follows from 

(61) 


Lemma 24 (ii). By Lemma 24 (i), 

dim (ker(S„,+,(C)^)f|5-+'(M^’"^')) +rs„^,(c)- = C\ 


m+Z 


In B.l and B.2 we prove that the summands in ( |61[ ) are bounded as 

dim (ker(S„+i(C)^) f| 




K-l 


, and 




respectively. Since - C™-^) = 

(iv) follow from (61). Statement (ii) follows from statement (i) and Lemma 25 below. 


C m —1 \ _ ^ 

R ! — ^ 


K+m+i-ij statements (i) and 

□ 
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B.l A lower bound on dim (keY{Sm+i{C)'^) f] 


In this subsection we prove the following result. 


Lemma 25. Let C G kc=K,m = R — K + 2,l>0, letF satisfy (301-(311, and let be defined 

by (36). Then 

(i) The matrix has full column rank, that is rp{m+i) = ; 


(ii) ker(S™+i(C)^)nS''"+'(K^"‘’'') D range(F(™+')). 

In particular, dim |^ker(Sm+z(C)^) p] > C^~^. 

Proof. Statement (i) was proved in [9l Proposition 1.10]. 

Let f be a column of the matrix F. Then f is a column of It is clear that G 

To prove (ii) we need to show that = 0. By Definition]^ the (ri,...,rm+/)th entry of the 

vector Sm is 


y (m + /)! 


E 


(si,...,s„+i)G 
-Pf-l -r^ + l} 





Since the vector f is orthogonal to exactly K — 1 columns of C, the fact that at least m indices of ri,..., rm+i 
are distinct, and m-\-l>R—K + 2 + l>R — {K — 1) it follows that (c^^f) • • • = 0, which completes 

the proof of (ii). □ 


B.2 A lower bound on the rank of the matrix S^+i (C) 

We need some additional notation. Let m > 2, I > 0, p > 0 and m<m + l—p<R. With an {m + l— p)-tuple 

(ii,..., im+i-p) such that 1 < zi < • • • < im+i-p < R 

we associate the set 


■— {(d) ■ • ■) im+i—p-! Iqi > • ■ • > dp) ' 1^9i — + Z p}. 

In other words, the set consists of the {m + Z)-tuples that are obtained by merging the {m + l — p)- 

tuple (ii,..., im+i-p) with p-combinations with repetitions of the set {ii,..., i 2 +i-p}. It is clear that for a fixed 
p there exist sets and each set contains (m +/(-tuples. 

Let E be the union of all sets Then the set E contains exactly 

C O I I I /^m ^m — 1 

l + l^R -h • " -h 


(to -I- /(-tuples (we follow the convention that ^ := 0 il m + I — p > R). Since, by construction, each 

(TO-|-/)-tuple of contains exactly m + l—p > to distinct elements, it follows that each (TO-|-/)-tuple of 

E contains at least to distinct elements. Let S^_|_;(C) denote the iL"* x matrix with columns 

(24), where (ri,... ,rm+i) € E. Then S^^i{C) is a submatrix of Sm+iiC). We have the following lemma. 


25 








Lemma 26. Let C G fcc = and m = R — K + 2. Then the matrix S^^;(C) has full column rank. 

In particular, rs„+,(c) > C'fl+/+i - . 

Proof. Suppose that there exists f G ^^r+i+i~^r such that S^_,_j(C)f = 0 . We show that f = 0 . We assume 
that the entries of f are indexed by (m + Z)-tuples (ri,... ,rm+i) G E, that is, in S^_,_;(C)f = 0 the column of 
S^+j(C) associated with the (m + 0-tuple (ii, ... ,im+i-p,iqi, ■ ■ ■ OgJ is multiplied by ■ 

To show that all entries are zero we proceed by induction on p = Z, ^ — 1,... ,max(0,m + 

I — R): in the pth step we assume that the identities 


=0, where 

1 < ii < • • • < im+i-p < R, l<<?i<---<<7p<2 + Z- p 

hold for p = Z,Z — 1 and prove that the identities hold for p = p. 

(i) Induction hypothesis: p = 1. We show that 


= 0 for 1 <ii < ■ ■ ■ < im < R, 1 < gi < • • • < ® < 2. 

We give the proof for the case L = = m, the other cases follow similarly. Thus, we show that 

. /=0f„rl £2. 

Since kc = K, the square matrix C := [ci C 2 Cm+i ■ • • c^j] is nonsingular. Let Ui and U 2 denote the first 
and the second column of C“^, respectively. Then 



1 0 
0 1 


0 ... 0 

0 ... 0 


( 62 ) 


Let X = tiUi +t2U2. Then the vector := x® ■ • -(gjx is orthogonal to the columns of the matrix S®_|_;(C) 

indexed by the (m + Z)-tuples (ri,..., rm+i) & E \ Ei . m- Indeed, if {ri,..., rm+i} \ { 1 , ■ • ■, m} 0 , then by 

pdl) and 


Hence 


^{m+l)T 

(to + ly. 


E 


‘'Sm + l 




X c 


■r,n+i. 


= 0 . 


( 63 ) 


(Si,.. 


rr+i)e 


0 = x(™+')^S^+;(C)f = (x^c^J • • • = 

'y ' (x Cl) ■ ■ ■ (x Cm)(x Cgj)---(x Cqj = ^ 04 ^ 

l<qi<--<qi<2 

(x^Ci) • • • (x^c„) Y (x'^C,J---(x^Cq,)/i„..,„,g, 

l<qi<--<qi<2 


Since fcc = K, at most one of the vectors Ui and U2 can be orthogonal to any of the vectors C3,..., Cm- Hence, 
(x^Ci) • • • (x^Cm) = tlt2(tiufc3 + t2U^C3) • • • (tiufc^ + t2U^C™) ^ 0 
for generic ti,t 2 G K. Hence, by ( [^ , 

Y (x^CqJ •••(x'^c,J/i,..=0 (65) 

l<qi<-<qi<2 
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for generic ti,t 2 € K. By construction of x, the I + 1 products 

(x^CqJ • • • (x^Cq,), 1 < qi < ■ ■ ■ < qi <2, coincide with the monomials Thus, identity (65) 


expresses the fact that a polynomial in ti and t 2 with coefficients vanishes for generic ti,t 2 £ ffi. 

It is well known that this is possible only if the polynomial is identically zero, yielding that fi,...^m,,qi,...,qi = 0 
for 1 < Qi < • • • < g/ < 2. 

(ii) Inductive step. We show that 


hi,. 




= 0 for 


1 < *1 < • • ■ < *r 


l + l — p ^ R-; 


l<qi< 


■ <qp <2 + 1—p 


or, equivalently, that = 0 for 

(*1) ■ ■ • ; im+l—pi iqi ) • • ■ ) iqp ) S 


u 


E.. 


—p ’ 




We give the proof for the case = 1,..., im+i-p = m + I — p, the other cases follow similarly. Thus, we show 
that fi,...,m+i-p,qi,...,qp = 0 for 1< qi<---<qp <2 + 1— p. The derivation is very similar to that of the 
induction hypothesis. 

Since fcc = K, the K x K matrix C := [ci ... C 2 +i-p c^+i-p+i ■ ■ ■ is nonsingular. Let Ui,..., U 2 +;_p 
denote the first 2 + I — p columns of C“^. Then 


r T 1 


T 

0 .. 

. 0 

0 . 

. 0‘ 


C = 

0 

1 .. 

. 0 

0 . 

. 0 

.^I+l-p. 


_0 

0 .. 

. 1 

0 . 

• 0_ 


Let X — tiUi + ■ • • + ^ 2 +/—pU 2 +;_p. Let 

Epcp := [J U 


( 66 ) 


and let the sets Ep^p and Ep=p be defined similarly. Then E = Ep^p U Ep^p U Ep=p. Then, by (24), (63), 
and (66), the vector := x ® ® x is orthogonal to the columns of the matrix S^_|_;(C) indexed by the 

(to + l)-tuples 

{ri,..., ^ m-\-l ) S Ep^p U {Ep—p \ ifi,.. .,m-\-l—p'} ■ 


Hence, similarly to (64) we obtain 


0 = x('"+')^S^+;(C)f = (x^c^J • • • = 

Y (x^C^J • ■ ■ (x^C^„+,)/ri.....r„,+,- 

( t *! 'j^Ep^pUEi ^ _ ^■mJj-l — p 
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Since, by the induction assumption, = 0 for (ri,... ,rm+i) S Ep^p, we have 


'y ' Cl) • • • (x Cm+i-p)(x Cgj)---(x = 

l<9l<'<9p<2+i—p 

(x Cl) • • • (x Cm+l—p) ^ ^ (x Cgj)---(x Cq^)fi^_^^rn+l—p,qi,...,qp- 

1<9i<---<9p< 2+/—p 


(67) 


Since fcc = Tf, at most 1 + Z — p of the vectors Ui,..., U 2 +i_p can be orthogonal to any of the vectors 
^s+i—pj ■ ■ • 7 ^m-\-i—p- Hence, 


(x^Ci) • • • (x^Cm+l-p) = 

tl • • • t2+l-p[tiU^ C3_|_/_p + • • • + t2+i-pU2+/_pC3_|_/_p) • • • 
(tlUi CjYi^i—p + * * ■ + ^2+/—pU2+/—pCm+/—p) 7^ 0 


for generic ti,..., t 2 +i-p G 


Hence, by (67), 


^ ' (x Cqj^)---(x Cqp)/i^,.._m+i-p,gi,...,ijp — 0 

l<9l<'<<?p<2+i—p 


( 68 ) 


(69) 


for generic ti,..., t 2 +z-p G K. By construction of x, the products (x^c^J • • • (x^Cq^), I < qi <■■■< Qp < 
2 + I — p, coincide with the monomials • • • t'^^i‘Sp}ai+--+a 2 +i-p=p- Thus, identity (69) expresses the fact 
that a polynomial in ti,..., t 2 +z-p with coefficients /i,...,m+z-p,iji,..., 9 p vanishes for generic fi,..., t 2 +z-p G It 
is well known that this is possible only if the polynomial is identically zero, yielding that fi,...,m+i-p,qi,...,qp = 0 
for 1 < qi < • • • < gp < 2 + 1 — p. □ 
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